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ABSTRACT 

We analyze the evolution of H II regions around the seven known SDSS quasars at z > 6. The comparison 
between observed and model radii of the H II regions generated by these quasars individually, suggests 
that the surrounding intergalactic hydrogen is significantly neutral. When all constraints are combined, 
the existing quasar sample implies a volume averaged neutral fraction that is larger than 10% at z > 6. 
This limited sample permits a preliminary analysis of the correlations between the quasar parameters, 
the sizes of their H II regions, and the associated constraints on the neutral hydrogen fraction. We find 
no evidence in these correlations to contradict the interpretation that the red side of the Gunn-Pcterson 
trough corresponds to the boundary between an H II region and a partially neutral IGM. 

Subject headings: cosmology: theory - galaxies: formation 



1. INTRODUCTION 

The Gunn-Peterson (1965, hereafter GP) troughs in the 
spectra of the most distant quasars at redshifts z ~ 6.3— 
6.4 (Fan et al. 2004) hint that the reionization of cosmic 
hydrogen was completed only a billion years after the big 
bang. Unfortunately, the troughs only set a lower limit of 
~ 10~ 3 on the volume averaged neutral fraction (Fan et 
al. 2001; White et al. 2003). This leaves open the ques- 
tion of whether reionization was completed only at z ~ 6 
so that the GP trough is due to a significantly neutral 
IGM, or whether it is only a residual neutral fraction left 
over from reionization at an earlier epoch that is respon- 
sible for the observed Lya absorption. This question has 
become particularly acute following evidence presented by 
the WMAP satellite in favor of early reionization (Kogut et 
al. 2003; Spergel et al. 2003). Indeed, if reionization ended 
only at z ~ 6, then more exotic "double reionization" sce- 
narios involving a massive early population of stars may 
be required (Wyithe & Loeb 2003a; Cen 2003). 

The H II regions generated by luminous quasars (Madau 
& Rees 2000; Cen & Haiman 2000) expand faster into a 
medium with a lower neutral fraction. Given information 
on the age of the quasar, we can then use the observed 
size of the H II regions to estimate the neutral fraction. 
In an earlier paper we have shown that the size of the 
H II regions around the two highest redshift quasars im- 
proves the lower limit obtained from the optical depth in 
the GP trough by two orders of magnitude [Wyithe & 
Loeb (2004a); see also supporting evidence by Messinger 
& Haiman (2004) who analyze SDSS 1030+0524 and find 
that the shape of the absorption spectrum near the edge of 
the H II region requires a large neutral fraction]. A large 
neutral fraction at z > 6 is in contradiction with the in- 
ference of a highly ionized IGM at z ~ 6.5 which has been 
made following the discovery of Lya emitting galaxies (e.g. 
Rhoads et al. 2004). However, the latter inference results 
from the assertion that the red-damping wing of resonant 
Lya absorption in the surrounding IGM would lead to a 
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strong suppression of the Lya line emitted by these galax- 
ies. This naive interpretation fails to account for clus- 
tering of sources which can generate large H II regions 
around the Lya emitters in an otherwise neutral IGM, 
and therefore allow transmission of the Lya line (Gnedin 
& Prada 2004; Furlanetto, Hernquist & Zaldariaga 2004; 
Wyithe & Loeb 2004c). 

In this paper, we supplement an earlier analysis of 
the H II regions around SDSS 1148+5251 and SDSS 
1030+0524 (Wyithe & Loeb 2004a) with new data on six 
additional quasars at z > 6. Our improved analysis incor- 
porates updated redshift data as well as the uncertainties 
in sizes of the H II regions, and in the ionizing luminosi- 
ties of quasars. In addition to the constraints derived from 
individual quasars, we compute combined constraints on 
the neutral fraction using all quasars. We also search for 
correlations between the quasar parameters, the sizes of 
their H II regions, and the associated constraints on the 
neutral fraction. 

The paper is organized as follows. In § 2 we summarize 
the observed properties of the sample, discuss the possibil- 
ities that a short Lya mean-free-path or dense absorbing 
cloud could mimic signature of an H II region, and con- 
sider the correlations between different observables. In § 3 
we describe our model for the evolution of an H II region 
surrounding a luminous z > 6 quasar. The limits that can 
be placed on the average neutral fraction of hydrogen in 
the intergalactic medium (IGM) around individual quasars 
are then described in § 4. Finally, we summarize our con- 
clusions in § 5. Throughout the paper we adopt the set of 
cosmological parameters determined by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP, Spergel et al. 2003), 
namely mass density parameters of VL m — 0.27 in matter, 
fib ~ 0.044 in baryons, fl\ = 0.73 in a cosmological con- 
stant, and a Hubble constant of Hq — 71 kms -1 Mpc -1 . 
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2. PROPERTIES OF THE KNOWN Z > 6 QUASARS 

The Sloan Digital Sky Survey (SDSS) has discovered 
seven quasars at z > 6 (Fan et al. 2001;2003;2004). These 
quasars are listed in Table 1 . Column 2 lists the UV mag- 
nitudes. Column 3 lists the source redshifts based on op- 
tically observed broad lines, such as Lya and high ion- 
ization lines such as N V, C IV, Si IV. Column 4 lists 
redshifts based on lower ionization lines observed in the 
near-IR, dominated by Mg II but including the Fe II com- 
plex (Maiolino et al. 2004; Frcudling et al. 2003; Iwamuro 
et al. 2004; Willott, McClure, & Jarvis 2003). Column 5 
lists the CO redshift for J1148+5251, and column 6 lists 
the redshift for the on-set of the GP absorption (see be- 
low). Column 7 lists the difference between the host galaxy 
redshift (see below) and the GP redshift, while column 8 
gives the corresponding size (i?) in physical Mpc of the 
H II region. 

In this paper we associate the red edge of the GP troughs 
in the spectra of z > 6 quasars with the presence of an 
H II region embedded in an IGM with neutral fraction 
Xhi. This assumption will be discussed further in §2.1. 
The key parameters in the analysis presented herein are 
therefore the host galaxy redshift and the redshift for the 
onset of GP absorption. Considering the host galaxy red- 
shift, the most accurate redshift estimate comes from the 
CO line emission from the host galaxy of J1148+5251, for 
which the uncertainty is Az — ±0.001 (Walter et al. 2003; 
Bertoldi et al. 2003). Note that accurate astrometry has 
shown that the CO emission is coincident with the optical 
QSO to within 0.1", providing strong evidence that the 
CO emission is from the QSO host galaxy (Walter et al. 
2004). Four other sources in Table 1 have reasonably ac- 
curate Mg II redshifts. Although Mg II is a broad line (few 
thousand km s _1 ), it is well documented to provide a rel- 
atively accurate estimate (within a few hundred km s _1 , 
or Az — ±0.01) of the host galaxy redshift in low redshift 
quasars (Tytler & Fan 1992; Richards et al. 2002). This 
expectation is supported by the good agreement between 
the CO and Mg II redshifts for J1148+5251. Two of the 
sources, J1623+3112 and J1602+4228, have only optical 
redshifts, corresponding to Lya and high ionization broad 
lines. The accuracy of the Lya emission line redshift is 
poor, due to strong associated absorption. Similarly, the 
high ionization lines are well known to be systematically 
blucshifted relative to their host galaxies, with a mean off- 
set of 824 ± 510km s _1 , as seen in a large sample of SDSS 
quasars (Richards et al. 2002). Hence, we assume a large 
error of Az = ±0.03 (^1200 kms -1 ) for these lines. In 
our analysis we rely on the CO, Mg II, and optical line 
redshifts, in descending order of preference. 

We estimate the GP redshift from the published spec- 
tra, as the redshift at which the emission extending blue- 
ward of the Lya (or Ly/3) peak becomes comparable to 
the noise. For one quasar, J1030+0524, the on-set of 
the Gunn-Peterson effect is well determined via sensitive 
Keck spectra of the Ly/3 region of the spectrum (White 
et al. 2003). A second source, J1148±5251, has a sensi- 
tive Keck spectrum in the Lya region (White et al. 2003), 
from which the on-set of GP absorption can be estimated 
fairly accurately. However, using Lya has the problem 
that the damping wing of a significantly neutral IGM just 
outside the Stromgren sphere can extend into the GP re- 



gion of the spectrum, thereby decreasing the apparent size 
of the sphere (Mesinger & Haiman 2004). For the one 
well documented case, J1030±0524, this effect amounts to 
Az = 0.02 (White et al. 2003), and we adopt this as the 
minimum error for the GP redshift based only on Lya. For 
most of the sources in Table 1 the on-set of the GP effect 
must be estimated based on the discovery spectra from 
3.5m telescopes, which acquire only a moderate signal-to- 
noise ratio. In these cases we again adopt a conservative 
GP redshift error of Az — 0.03. 

In addition to uncovering the onset of a GP trough, the 
deep spectrum of SDSS 1148+5251 shows a peak of trans- 
mitted flux which is observed both in the Lya and Ly/3 
troughs, indicating the presence of a transparent window 
in the IGM (White et al. 2003). More puzzling is the ob- 
servation of transmission peaks in the Ly/3 trough that 
do not have counterparts in the Lya trough. White et 
al. (2004) interpret this as an indication of the presence of 
two overlapping windows of transparency; one in the Ly/3 
forest at z ~ 6, and one in the Lya forest at z ~ 5. Al- 
though this scenario is unlikely a-priori, White et al. (2004) 
point out that HST imaging did not reveal a small cluster 
of galaxies at z ~ 5 as had originally been hypothesized. 
It is still possible that a single galaxy, whose angular po- 
sition is coincident with the quasar could generate both 
the H II region required for transparency, as well as the 
flux in the transmitted peak. Such a galaxy could also ex- 
plain the presence of low level continuum in the Lya and 
Ly/3 troughs. However as pointed out by Oh & Furlan- 
etto (2004), transmission of continuum is also detected in 
the Ly7 trough, which argues against a foreground galac- 
tic contribution which would be subject to a Lyman break. 
Instead Oh & Furlanetto (2004) argue that transmission 
of continuum in the Lya, Ly/3 and Ly7 troughs points to 
a highly ionized IGM at z ~ 6.2 along the line of sight to 
this quasar. 

Estimates of the neutral fraction based on transmission 
in the Lya and Ly/3 troughs are made at a redshift in the 
center of the trough, while estimates based on the size of 
the H II region are made at the redshift of the quasar. For 
the purpose of the present work we assume that the red 
edge of the Lya and Ly/3 troughs in SDSS 1148±5251 mark 
the boundary of an H II region (this point is discussed fur- 
ther in §2.1). Along any line of sight, the redshift of Lya 
transmission will be somewhat lower than the redshift at 
which the IGM became reionized. This latter redshift is 
expected to have a scatter among different lines-of-sight of 
at least 0.15 (Wyithe & Loeb 2004d). It would therefore 
not be surprising if the IGM were reionized at z ~ 6.2 
along the line of sight to SDSS 1148+5251, but only at 
z <~ 6 along the lines of sight to other quasars. 

While better optical and near-IR spectra, as well as 
more CO host galaxy redshifts, would improve the analy- 
sis presented herein (and such programs arc in progress), 
we already find that interesting conclusions can be drawn 
even for the very conservative errors adopted in Table 1. 
Given that other uncertainties in the problem (e.g. the 
ionizing luminosities and lifetimes of the quasar) are of 
comparable importance, the precise determination of red- 
shifts will not affect significantly the statistical conclusions 
reached in this study. 
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Table 1 
Properties of z > 6 quasars 



Source 


Mi450 


ZLya,CIV 


ZMgll 


zco 


ZGP 


Zhost — ZQP 


R±o R 


J^(R) 


J1148+5251 


-27.82° 


6.37±0.03 c 


6.403 ± 0.005 9 


6.419 ±0.001" 


6.325 ± 0.02 c 


0.09 ± 0.02 


4.89±1.09 


1.0 


J1030+0524 


-27.15 d 


6.28±0.03 c 


6.311 ±0.005 9 




6.178 ±0.005 c 


0.133 ±0.007 


7.50±0.39 


0.2 


J1048+4637 


-27.55 a 


6.23 ±0.03° 


6.203 ± 0.005 9 




6.16 ±0.03° 


0.04 ±0.03 


2.34±1.76 


3.4 


J1623+3112 


-26.67 e 


6.22 ± 0.03 e 






6.16±0.03 e 


0.08 ±0.04 


4.65±1.74 


0.4 


J1602+4228 


-26.82 e 


6.07±0.03 e 






5.95±0.03 e 


0.12 ±0.04 


7.36±2.45 


0.2 


J1630+4012 


-26.11° 


6.05 ±0.03° 


6.065 ± 0.005 s 




5.98 ±0.03° 


0.085 ±0.03 


5.22±1.84 


0.2 


J1306+0356 


-27.19 d 


5.99 ± 0.03 d 


5.99±0.02 / 




5.90 ± 0.03 d 


0.09 ±0.04 


5.68±2.52 


0.4 



"Fan et al. (2003); "Walter et al. (2003); "Bertoldi et al. (2003); c White et al. (2003) 
d Fan et al. (2001); e Fan ct al. (2004); / Maiolino et al. (2003); 9 Iwamuro et al. (2004) 



2.1. Could the Red Boundary of the GP Trough 
Correspond to the Lya Mean- Free- Path in a Highly 
Ionized IGM Rather than an H II Region? 

Observations of the Lya forest show that the number of 
absorbing systems increases towards high redshift. This 
results in a decrease of the mean-free-path for ionizing 
photons (Fan et al. 2001; Miralda-Escude 2003) and an 
increase in the optical depth for Lya transmission (i.e. 
the GP optical depth) with increasing redshift. Fan et 
al. (2001) find that the mean-free-path for ionizing pho- 
tons was ~ 7 physical Mpc at z ~ 5.5 (comparable to the 
size of the quasar H II regions in our discussion), but had 
declined below an upper limit of ~ 1 physical Mpc within 
the GP trough along the line of sight to SDSS 1030±0524. 
In this context, the IGM transmission just blueward of the 
Lya resonance of the z > 6 quasars extends out to ~ 5 
physical Mpc away from them, and must be caused by 
their ionizing radiation. A natural interpretation is that 
these quasars ionize the surrounding IGM to a level well 
in excess of that expected from an extrapolation the ioniz- 
ing background from lower redshifts. Indeed, the ionizing 
flux from the quasar SDSS 1030+0524 should exceed the 
estimated ionizing background at z ~ 5.5 out to a distance 
of ~ 10 physical Mpc (at z ~ 6.2). These considerations 
imply that the mean-free-path of ionizing photons should 
be at least as large as 5-10 physical Mpc around the z > 6 
quasars. 

Now suppose the IGM were already highly ionized at 
z > 6. One may wonder whether the red edge of the GP 
trough does not actually correspond to the boundary of an 
H II region, but rather to the radius where Lya photons 
near the quasar are absorbed by residual H I in an other- 
wise highly ionized IGM. We have argued that in a highly 
ionized IGM one could expect the ionized fraction of the 
IGM within - 10 Mpc of the quasar SDSS 1030±0524 to 
exceed that at z ~ 5.5. This ionization state would be con- 
trary to the observation of a GP trough, which points to 
a rapid change in the ionization state of the IGM at z ~ 6 
(Fan et al. 2001). Nevertheless, we would like to critically 
examine the possibility of a highly ionized IGM and Lya 
absorption from residual neutral hydrogen. To this end 
we find whether the ionizing intensity of the quasar at the 
edge of the HII region is in excess of that required to lower 
the GP optical depth below the observed limit. 

Table 2 summarizes the absorber redshift z a b s (col- 
umn 1), effective GP optical depth r c ff(z a b s ) at z a b s (col- 



umn 2), and background ionizing intensity </2i( z abs) at 
z a bs in units of 10 -21 erg s _1 Hz -1 cm~ 2 sr _1 (column 3) 
as inferred from spectra of high redshift quasars (Fan et 
al. 2001). At a fixed level of dumpiness in the IGM, the ef- 
fective GP optical depth r ff (zqp) due to neutral hydrogen 
in ionization equilibrium with the quasars ionizing radia- 
tion field at zqp follows 

n- (y \ t ( 9 \(1±I2L ^2l(£abs) m 

r off (z GP ) ~ r cff (z abs ) j (1) 

(see Eq. 6 in Barkana & Loeb 2004), where ^'(fl) 
is the quasars ionizing intensity at a distance R 
from the quasar (column 8 of Table 1) in units of 
10 -21 erg s _1 Hz -1 cm~ 2 sr _1 . We use a specific ionizing 
luminosity of 10 18 05 ergs _1 Hz -1 per solar B-band lumi- 
nosity (Schirber & Bullock 2002). These values are listed 
in column 9 of table 1. From equation (1) we evaluate 
the effective optical depth T c ff(zGp) a t redshift zqp for 
each quasar, extrapolating from absorbers at z a b s ~ 5 and 
z a bs ~ 5.5 respectively. The values are listed in columns 
4-10 of table 2 for the z > 6 quasars considered in this 
paper. 

The extrapolation of optical depth to redshifts beyond 6 
using equation (1) assumes a clumping factor that does not 
evolve in time. However clumping in the IGM is expected 
to increase towards low redshift, as structure grows in the 
IGM. For this reason the estimates of the effective GP opti- 
cal depths surrounding the high redshift quasars are larger 
if extrapolated from absorbers at lower redshift. The most 
appropriate estimates of optical depth due to neutral hy- 
drogen in ionization equilibrium with the quasar flux at a 
distance R are therefore given by the second row in Ta- 
ble 2. These values lie in the range of r e e = 0.1 — 1.9, which 
would allow transmission in the observed GP troughs out 
to distances that are well in excess of ~ 5Mpc away from 
the quasar. If the quasars were emitting into an already 
ionized IGM, then the spectra would show the usual prox- 
imity effect, with a Lya forest persisting at redshifts where 
the GP troughs are actually observed. We therefore infer 
that the onset of the GP trough is caused by diffuse neu- 
tral IGM rather than the mean free path of Lya photons 
near the quasar. 

2.2. Could the Lya Transmission be Terminated by a 
Dense Clump, Rather than the Neutral IGM? 

Our analysis of limits on the neutral fraction from the 
size of H II regions relies on the inference that the red end 



4 



Table 2 

Lyq Absorption Properties Near z > 6 Quasars 



Z&hs 


Tcff (.Zabs) 


•^21 (Zabs) 


Extrapolated Lya Optical Depth Surrounding z > 6 Quasars 
J1148 J1030 J1048 J1623 J1602 J1630 J1306 


5±0.2 
5.5±0.2 


2.0 ±0.3 
2.5 ±0.3 


0.1±0.02 
0.1±0.02 


0.5 1.9 0.1 1.1 2.2 2.1 0.9 
0.4 1.7 0.1 1.0 1.9 1.9 0.8 



of the GP trough can be used to estimate the size of H II 
regions. However the spectra of the z > 6 quasars show 
the presence of clouds with a substantial column density. 
If there had been a cloud of sufficiently high column den- 
sity to produce an absorption wing over a sufficiently large 
wavelength range, the H II region would have appeared 
smaller than it actually is, resulting in lower limits that 
are systematically high. We believe that this is unlikely 
for the following reasons. 

We start from an observational perspective. The hy- 
drogen column density required to produce a Lya optical 
depth greater than 26 [observed limits for the best stud- 
ied example (White et al. 2003)] across a significant frac- 
tion, e.g. 10%, of the 100A within the observed ionized 
region is N = 10 21 cm~ 2 at this rcdshift. This column 
density corresponds to damped Lya absorbers, which are 
extremely rare (Storrie-Lombardi & Wolfe 2000), ~ 0.1 
per unit redshift between z ~ 2 and z ~ 4. Thus it is very 
unlikely that a damped Lya absorber would terminate the 
GP trough in any one of the 7 quasars, as we are dealing 
with path lengths of only ~ 5 Mpc (or a redshift interval 
of&z~0.1). 

This small probability can also be substantiated by a 
simple theoretical argument. Our model adopts the dis- 
tribution of gas clump overdensitics from the numerical 
simulations of Miralda-Escude et al. (2000) to derive the 
critical overdensity, A cr ;t, such that the typical line of sight 
to the edge of the H II region at ~ 5 Mpc would not con- 
tain any clumps with overdensity larger than A cr it, while 
all clumps with overdensity below A cr i t are ionized. We 
find that A cr it = 20. Barkana & Loeb (2002) have proven 
a simple theorem stating that: "The fraction of the line- 
of-sight covered by gas at a given overdensity is equal to 
the volume filling fraction of gas at that overdensity". We 
may therefore find the total column density in clouds along 
the line-of-sight with overdensities greater than A cr i t as 
follows. The fraction of mass contained within overdensi- 
ties greater than A cr it (~ 0.07) is multiplied by the mean 
column density of neutral gas along a path through the 
neutral IGM of length equal to the size of the H II region 
(1.4 x 10 21 cm~ 2 ). We then find the total column den- 
sity due to clouds along the line-of-sight with overdensities 
greater than A crit to be 0.07 x 1.4 x 10 21 = 10 20 cm" 2 . This 
total column density falls short by an order of magnitude 
relative to the value required to produce a wide wing of 
absorption that would affect significantly the size estimate 
of the H II region. 

2.3. Distributions of properties among z > 6 quasars 

We now examine the correlation between the properties 
of the z > 6 quasars and the sizes of their H II regions. 
The left hand panel of Figure 1 shows the correlation be- 



tween the measured absolute magnitude at a rest-frame 
wavelength of 1450A, M1450, and the inferred radius of 
the H II region, R. The line corresponds to the naive ex- 
pectation that R oc L 1 / 3 (normalized to the quasar with 
the largest value of RL~ X / 3 ). Scatter would be expected 
around this line for physical reasons (e.g. owing to differ- 
ent quasar ages) in addition to measurement errors. How- 
ever, comparison of this line with data points constitutes 
an important systematic check that the bubbles do not 
show any unexpected trend. 

Since the volume within the H II region grows approxi- 
mately linearly with time, the observed bubble sizes should 
follow a characteristic distribution. The central panel of 
Figure 1 shows a cumulative probability histogram for 
T] = RL^ 1 / 3 (normalized to the quasar with the largest 
value of RL~ X / 3 so that < r] < 1). For a random sample 
of quasar ages, this distribution would be expected to fol- 
low P(< rj) = rj 3 (solid line). The expected and measured 
curves are fully consistent in a Kolmogorov-Smirnov (KS) 
test with a probability ffcs ~ 0.7. On the other hand if 
the observed size of the Lya transmission region had been 
set by a dense absorbing cloud rather than by the neutral 
IGM, then we would expect a random distribution of bub- 
ble sizes R. Since R and L would be uncorrelated in this 
case, the distribution of RL^ 1 ^ 3 would also be random. 
However we find that a random distribution (dashed line) 
is also consistent in a KS test, with Pks ~ 0.2. 

We note one possible anomaly in the rcdshifts listed in 
Table 1. The object SDSS 1048±4637 is the only source 
with a value of ZMgii < ZLya.civ, which is not expected 
given the distribution of blucshifts for high ionization lines 
(Richards et al. 2002). SDSS 1048±4637 is also the only 
object with a value of ZMgii — ^gpt, and hence a size 
for the H II region that is consistent with zero (open 
circle in the left hand panel of Figure 1). One would 
need an observationally motivated reason to eliminate this 
quasar from the sample; nevertheless, it is interesting to 
repeat the above comparison of distributions in the ab- 
sence of SDSS 1048±4637 (right hand panel of Figure 1). 
In this case we see that while the distribution is consistent 
with volumes of H II regions that grow linearly with time 
(-Fks ~ 0.9), a random distribution is now only consis- 
tent with the data at the 7% level. In the future, using 
a larger sample of quasars (perhaps from the completed 
SDSS) we may be able to reject the random distribution. 
For now, we will return to the possibility that dense ab- 
sorbing clouds might mimic the signature of H II regions 
in a neutral IGM in section 2.2. 

Another possible source of systematic uncertainty arises 
from gravitational lensing. The large magnification bias 
for quasars on the bright end of the quasar luminosity 
function implies that strong gravitational lensing may be 
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Fig. 1. — Correlations between the properties of the z > 6 quasars and their H II regions. Left panel: correlation between quasar 
absolute magnitude (at rest-frame 1450A) M1450 and inferred radius of the H II region around it R. The open circle shows the point 
for SDSS 1048+4637. The line corresponds to the naive expectation that R oc L 1 / 3 . Central panel: cumulative probability histogram for 
r; = RL~ 1 / 3 /(RL~ 1 / 3 ) max . The solid line corresponds to P(< rj) = rp , and the dashed line to a random distribution. Right panel: same as 
for the central panel, but with SDSS 1048+4637 omitted. 



one or even two orders of magnitude more common in 
the SDSS z ~ 6 quasars than in lower redshift samples 
(Wyithe & Loeb 2002a;2002b; Commerford, Haiman & 
Schaye 2002). Undetected lensing leads to an overesti- 
mate of the quasar luminosity. In the present analysis 
this leads to an overestimate of the neutral fraction. Cur- 
rently there is high resolution imaging for four of the z > 6 
quasars. As part of an HST snapshot survey of high red- 
shift quasars, SDSS 1030+0524 and SDSS 1306+0356 have 
been imaged (Richards et al. 2003). In addition Keck K- 
band images have been obtained for SDSS 1048+4637 and 
SDSS 1148+5251 (Fan et al. 2003), while SDSS 1148+5251 
has also been imaged with HST (White et al. 2004). None 
of these images show evidence for strong lensing, and hence 
for significant magnification (however see Keeton, Kuhlun 
& Haiman 2004). As a result we do not consider lensing 
in our analysis. 

3. THE MODEL 

The model used in this paper to compute the evolution 
of quasar H II regions has been previously described in 
Wyithe & Loeb (2004a). However for completeness we 
present a summary of its main features below. 

Quasars are assumed to be powered during the hierar- 
chical growth of their host galaxy. Bright episodes are trig- 
gered when halos merge. Within our model, we generate 
many random realizations of the merger tree of the host 
halo at z = 6.5, assign super-massive black holes (SMBHs) 
to these halos and compute the time-dependent luminosity 
that is triggered during the merges. Each tree has N melge 
major mergers, which occur at times tj. We assume a re- 
lation between the black-hole and galaxy halo mass that 
preserves the correlation between the circular velocity of 
the halo and the black-hole mass it hosts in the form 



M halo = 1.5 x 1O 12 M 



Mi 



bh 



3/5 



1 



7 



-3/2 



(2) 



V 10 9 M Q / 

Inferred black-hole masses for the SDSS quasars of Mbh ~ 
1O 9 M imply halo masses of M halo ~ 1O 12 M (Wyithe 
& Loeb 2003). Since this relation is non-linear, there is 
a mass differential between the coalesced black-hole mass 
and the mass of the new halo. We define a parameter 
fit = /mass/??Edd, where / mass is the fraction of the mass 



differential that is accreted during the luminous phase and 
?7Edd is the fraction of the Eddington rate at which the 
mass is accreted. If a merger of two halos leads to coa- 
lescence of their black-holes, and the mass differential is 
accreted during a luminous phase over which the quasar 
shines near its Eddington limiting luminosity, then the 
quasar lifetime is 

~(M 1 +M 2 ) 5/3 ~ 



tit = fit x 4 x 10' 



0.1 



In 



5/3 



M. 



5/3 



years, (3) 



where M\ and Mi are the masses of the merging halos, and 
e (taken to be 0.1) is the radiative efficiency. Note that 
the lifetime increases for sub-Eddington accretion. The pa- 
rameter /u may be thought of as the fraction of our fiducial 
lifetime during which the quasars shine, leading to quasar 
lifetimes of t\ t ~ fit x 10 7 years. The model lifetime is 
therefore consistent with current estimates [10 6 — 10 s years; 
see Martini (2003) for a summary] for values of /it that are 
of order unity. We assume that each quasar episode j has 
an exponential light-curve 

Nj (t) = e(t-tj) No j cxp {-(t-tj) I tit } (4) 
beginning at tj and with a characteristic decay time of 
tn during which the SMBH shines at its Eddington lumi- 
nosity, L E = 1.4 x 10 47 erg s- 1 (M bh /lO 9 M ). Here 9 
is the Heaviside step function. The time dependent ion- 
izing luminosity within the merger tree is then computed 
from the sum of luminous episodes N(t) = y~*^™ TS " j\r^ 
Following Telfer et al. (2002) and White et al. (2003), 
we adopt an ionizing photon emission rate of Nqj = 
2 x 10 57 s -1 / (aEUV,j/ — 1-75), where «EUV,j is the spec- 
tral index in the EUV band during the jth merger. Telfer 
et al. (2002) find o:euv = —1-75 ± 0.75, and we assign a 
value of «EUV,i to each merger in the tree from a normal 
distribution of mean -1.75 and standard deviation 0.75. 

The evolution of the physical radius of the H II region, 
R, may then be found through integration of the differen- 
tial equation 



dR 
~dt 



F^Njt) - a B CF m x m (n*y (f R 3 ) 
F 7 N(t) + AttR 2 (1 + zY 1 cF m x m n% 



(5) 



where c is the speed of light, is the mean number 
density of protons at z = 0, cub = 2.6 x 10~ 13 cm 3 s -1 
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is the case-B recombination coefficient at the characteris- 
tic temperature of 10 4 K, and Ni on is the rate of ionizing 
photons crossing a shell at the radius of the H II region 
at time t. We use the distribution derived from numer- 
ical simulations for the over-densities A in gas clumps 
(Miralda-Escude, Haehnelt & Rees 1998), and calculate 
the mean free path for ionizing photons d(A c ) as a function 
of the critical overdensity A c (Miralda-Escude, Haehnelt 
& Rees 1998; Barkana & Loeb 2002). Following Barkana & 
Loeb (2002), we then find the value of A c at which a frac- 
tion exp [— R/d(A c )] — 50%, of the emitted photons do not 
encounter an overdensity larger than A c within the H II 
region. We also compute the mass fraction F m (~ 1) of gas 
within R that is at over-densities lower than A c . Finally, 
we calculate the clumping factor in the ionized regions, 
C(R) = (A 2 )/(A) 2 , where the angular brackets denote an 
average over all regions with A < A c . For xhi = 1 and 
dR/dt >C c, equation (5) reduces to its well-known form 
(e.g. Madau & Rees 2001). The emission rate of ionizing 
photons N in equation (5) is computed at t' = t — tdday to 
account for the finite light travel time between the source 
and the ionization front. 

Note that while equation (5) is expressed in terms of the 
radius of a spherical H II region, there is no implicit as- 
sumption about isotropy in the analysis presented in this 
paper. This is because both the volume of the HII region 
and the luminosity of the quasar are measured per unit 
solid angle along the line-of sight. The extrapolations to 
total volume and total luminosity are made by multiplying 
these quantities by 4ir purely for convenience of presenta- 
tion. 

In our calculation we require that the quasar ionize re- 
gions up to a sufficiently high overdensity, so as to allow 
the ionizing photon mean-free-path to exceed the radius 
of the H II region. The highest density regions (which are 
sufficiently compact to allow a long mean-free-path) may 
remain neutral. The value of xm in equation (5) refers to 
the neutral fraction in the low density regions. We there- 
fore interpret the neutral fraction xm in equation (5) as 
volume weighted rather than mass weighted fraction. 

In the hierarchical picture of structure formation, the 
appearance of the quasar and the surrounding galaxies will 
occur concurrently. The neutral fraction into which the 
quasar H II region expands therefore reflects the contri- 
bution to reionization due to stellar flux from the quasar 
host and surrounding galaxies (see Fig. 2 in Wyithe & 
Loeb 2004c for their relative significance). As the quasar 
and stellar ionizing fluxes are emitted at the same cosmic 
epoch, we implicitly assume that both are responsible for 
reionizing the low density regions, such that both have an 
ionizing photon mean free path that exceeds the radius of 
the H II region. In this our model differs markedly from 
the work described in Yu & Lu (2004), where quasar flux 
is assumed to be emitted into an IGM whose low density 
regions have already been ionized by stars at some prior 
epoch. In that work the clumping factor associated with 
quasar ionizing flux is evaluated above the density thresh- 
old corresponding to a mean-free-path that equals the size 
of the H II region, i.e. it is assumed that stars ionize the 
low density regions while the quasar ionizes only the high 
density regions. As a result, Yu & Lu (2004) infer a much 
higher clumping factor than we find here, and conclude in 



difference to this work that the quasar flux does not pro- 
vide a significant contribution to the growth of the H II 
regions. 

Several consequences of our simple model can be imme- 
diately confronted with observational data. Evidence from 
direct determination of the Mbh-velocity dispersion rela- 
tion (Shields et al. 2003), as well as the redshift evolution 
of the quasar correlation function (Wyithe & Loeb 2004b) , 
suggest that the assumed relation (equation 2) is indeed 
preserved out to a redshift of at least z ~ 3. Prelimi- 
nary support for an extrapolation of the relation out to 
z ~ 6 comes from the velocity shift of the Lya absorption 
feature due to the galactic virialization shock, whose am- 
plitude gauges the circular velocity and hence mass of the 
host halo for some high-redshift quasars (Barkana & Loeb 
2003). The presence of a massive galaxy is also implied 
by the molecular mass of > 10 10 M© in the host galaxy 
of SDSS1 148+5251 (Walter et al. 2003; 2004) and the ve- 
locity width of 280 km s _1 for its CO lines (Bertoldi et 
al. 2003). The measured CO velocity width corresponds 
to a dark halo mass of ~ 10 11 M Q at z = 6.4, which is 
insufficient to contain the inferred molecular gas mass. In- 
deed, 280km s -1 is much smaller than the ~ 500km s _1 we 
would expect to be associated with a SMBH of ~ 1O 9 M 
(Wyithe & Loeb 2003,2004b), which is thought to power 
the observed quasar. A possible explanation for this in- 
consistency is that the CO observations only sample the 
gravitational potential within a few kpc from the galaxy 
center while the galaxy halo has a much deeper potential 
well. To quantify this uncertainty we note that the num- 
ber of major mergers per logarithm of mass increases by 
a factor of ~ 2.4 if the host dark matter halo mass for the 
SDSS quasars is an order of magnitude smaller than as- 
sumed (i.e. 10 n M Q rather than 10 12 M Q ). Thus, adopting 
a smaller dark matter halo at the bottom of the merger 
tree would lead to a higher merger rate, a more frequent 
quasar activity, and hence the inference of a larger neutral 
fraction. 

Additional support for the assumed Mhaio— relation 
comes from attempts to model the luminosity function of 
quasars using the abundance of halos and the Mbh — o re- 
lation (Volenteri, Haardt & Madau 2002; Wyithe & Loeb 
2003). These models are equally successful at z ~ 6 as 
they are at lower redshifts z ~ 2-3 where data exists on 
the Mbh — cr relation. These physically motivated models 
for the luminosity function also suggest a quasar lifetime 
that is in the vicinity of 10 6 — 10 7 years at z ~ 6. Com- 
paring with tit (as calculated from equation 3) implies a 
value f\t > 0.1, or in other words that most of the SMBH 
mass at z > 6 was accreted during the luminous phase. 
Indeed having all the black-hole mass accreted during the 
luminous phase (/i t = 1) is consistent with the census of 
mass in local dormant SMBHs, compared with the mass 
accreted during luminous quasar phases throughout the 
history of the universe. These studies are most sensitive 
to conditions at z ~ 2 — 3, but find the majority of SMBH 
mass to have been accreted during luminous quasar phases 
near the Eddington limit, and with a radiative efficiency 
of e - 0.1 (Yu & Tremaine 2002). 

To estimate the ionizing luminosity of the z > 6 quasars 
we use the median spectrum of low redshift quasars derived 
by Telfer et al. (2002), scaled to the appropriate luminosity 
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of the quasar (M1450). The assumption of no evolution of 
the quasar spectrum over 90% of the age of the universe is 
supported by the observation that the median rest frame 
UV spectrum of the high redshift quasars is consistent 
with that at low redshift (Fan et al. 2004). In addition, 
the recently observed X-ray spectrum of SDSS 1306+0356 
implies an optical to X-ray spectral index that is consistent 
with radio-quiet quasars at lower redshift (Schwartz & Vi- 
rani 2004). These results imply little evolution in quasar 
spectra, and justify our use of the low redshift median 
quasar spectrum for the analysis of the z > 6 quasars. 

4. CONSTRAINTS ON THE NEUTRAL FRACTION 

Next we derive limits on the neutral fraction of hydro- 
gen surrounding the z > 6 quasars. Our basic method 
follows that outlined in Wyithe & Loeb (2004a). However 
we have augmented the approach to reflect the richer data 
set now available. 

We compute the conditional probability distributions 
^kLhi ^ or ^ e 0Dserve d radius of the H II region sur- 
rounding each quasar as a function of the neutral fraction 
of hydrogen, xm- Based on these distributions we find the 
likelihood for the neutral fraction from the observed radius 
around each quasar 

Li(xm) = [ dRiN(R?», a Ri ) ^(R = Ri) (6) 

J ^ XHI 

Here N(x,a x ) is a normal distribution of mean x and vari- 
ance a x . The values of mean and variance of the radii are 
given in Table 1. Assuming the neutral fraction to have 
the same value around all quasars in the sample, we also 
find the joint likelihood 

L(x m ) = Ui=i,7Li. (7) 
The relative likelihood for xm may be combined with 
its a-priori probability distribution dP pr i OI /dxm to yield 
cumulative a-posteriori probability distributions 



Pi(< x m . 



Jo 



, , APprior 

dx m Li(x m )- 



dx 



HI 



and 



P{< x m ) = N 



dP ■ 
j I ri I \ prior 



dx[ 



(8) 



(9) 



HI 



where N and TVj are normalizing constants such that 
Pi(< 1) = 1 and P(< 1) = 1 respectively. In this paper 
we adopt a logarithmic prior for XHI; /dxm oc 

1/xm for 10~ 3 < xm < 1, corresponding to the range 
allowed by observations of the GP trough (White et 
al. 2003). More stringent limits are found for xhi if a 
flat prior is assumed. 

In Figure 2 we plot the cumulative probability of 
xhi- The thin lines show the distributions for individual 
quasars, whereas the thick line shows the combined result. 
The left panel shows the case of f\ t = 0.1, correspond- 
ing to the situation where the quasar fiducial lifetime in 
equation (3) overestimates the true lifetime by a factor of 
10. The right hand panel shows the fiducial case, with 
/it = 1. We note that each of the seven quasars individ- 
ually yields a consistent result that xm is of order 0.1-1. 
The combined constraint shown by the thick line is there- 
fore not controlled by observations of just one or a couple 
of the quasar H II regions. The combined constraint for 



the fiducial model is xm > 0.1 (fit = 0.1) and xm > 0.7 
(/it = 1) at the 90% confidence level. Clearly /i t is a limit- 
ing systematic uncertainty and we return to its dependence 
below. 

In Figure 3 we examine the dependence of the individual 
quasar limits on the quasar redshift. Each quasar is repre- 
sented by 4 points (with SDSS 1048+4637 denoted by open 
circles), which from bottom up show the 1st, 10th, 90th 
and 99th percentiles of Pj(< xhi)- Lines join these points 
for each percentile to guide the eye. The left and right- 
hand panels illustrate the cases of /i t =0.1 and j\ t = 1, 
respectively. The apparent trend is that higher limits on 
xm are derived for higher redshift quasars. This trend is 
to be expected if the quasars are observed near the end of 
the reionization process. 

The variation of constraints with quasar redshift may 
be compared with expectations from theory. The redshift 
at which a region of a given size is reionized is propor- 
tional to the linear overdensity on that scale (Barkana 
& Loeb 2003). Along different lines of sight this red- 
shift has a Gaussian distribution with a variance given by 
the power spectrum of primordial fluctuations. Wyithe & 
Loeb (2004d) have shown that the combined constraints of 
finite light travel time and cosmic variance imply that the 
scatter in the redshift at which neutral IGM would be en- 
countered along a random line of sight is Az ~ 0.15. This 
scatter defines a minimum rate over which the IGM could 
become reionized. For comparison with the observed lim- 
its, the dashed line in Figure 3 therefore shows the cumu- 
lative distribution corresponding to a Gaussian with vari- 
ance 0.15 around a central redshift of 6.4 (which results in 
a neutral fraction of xm ~ 10~ 3 at z ~ 5.85 where limits 
exist based on the optical depth in the GP trough). The 
limits from individual quasars imply a somewhat slower 
evolution in the neutral fraction than the maximum rate 
described by the dashed curve, as expected. 

As already mentioned, the largest systematic uncer- 
tainty in the analysis described is the value of f\ t . While 
we expect /it ~ 1, it is instructive to compute the lim- 
its on the neutral fraction as a function of /i t . In Fig- 
ure 4 we show the 1st, 10th, 90th and 99th percentiles of 
P{< xm) from the bottom up. Smaller values of /it yield 
less stringent limits. However we find that Xm > 0.1 for 
all /it > 0.1, and x m > 0.01 for all / it > 0.025 at 90% con- 
fidence. These limits represent a significant improvement 
upon the limits for the volume averaged neutral fraction 
based on the optical depth in the GP trough, xm > 10~ 3 
(White et al. 2003). 

5. CONCLUSION 

In this paper we have extended an earlier analy- 
sis of the neutral fraction of hydrogen in the IGM 
around SDSS 1148+5251 and SDSS 1030+0524 (Wyithe 
& Loeb 2004a) to include all seven quasars now known at 
z > 6. We have used updated redshifts for the hosts of 
these quasars and incorporated uncertainties in the mea- 
sured size of their H II region. The small size of the H II 
regions implies that the typical neutral hydrogen fraction 
of the IGM beyond z ~ 6 is in the range 0.1-1. This re- 
sult also holds for the IGM surrounding each individual 
quasar when the H II regions are considered separately. 
The primary systematic uncertainty in the analysis is the 
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Fig. 2. — Cumulative probability for the neutral fraction of the IGM, xm. The thin lines show the distributions for individual quasars, 
while the thick line shows the combined result. Left: f\ t = 0.1. Right: f\ t = 1. 




Fig. 3. — Limits from individual quasars as a function of their redshift. Each quasar is represented by a point. From bottom up, the curves 
show the 1st, 10th, 90th and 99th percentiles of P(< 2?hi)- The left and right-hand panels show the cases of /i t = 0.1 and f\ t = 1, respectively, 
and limits derived from SDSS 1048+4637 arc denoted by an open circle. The dotted curves provide the integral of a normal distribution in 
redshift Jgdz'N (6.4, 0.15). 
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Fig. 4. — The neutral fraction as a function of /it. From the bottom up, the lines show the 1st, 10th, 90th and 99th percentiles of P(< ^hi) 
respectively. 
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quasar lifetime. However by combining the limits for the 
six quasars, we find that at the 99% level, the neutral frac- 
tion is larger than 0.08 for quasar lifetimes > 10 6 years, or 
larger than 0.6 for lifetimes > 10 7 years. These lifetimes 
correspond to only 10~ 3 and 10~ 2 Hubble times respec- 
tively at z ~ 6 and are favored by empirical constraints on 
the lifetime of lower-redshift quasars (e.g. Martini 2003). 
Larger duty-cycles lead to stronger limits. We find that 
the size distribution of H II regions is consistent with the 
expected distribution for observations that are made ran- 
domly in time. In addition, we find that constraints on 
the neutral fraction obtained from individual quasars are 
stronger at higher redshift. This is to be expected if the 



universe is nearing the end of reionization at z ~ 6 as 
the neutral fraction drops with time. A larger sample 
of quasar H II regions from the full SDSS, or those dis- 
covered by forthcoming redshifted 21cm surveys (Wyithe 
& Loeb 2004e) will allow more stringent checks of these 
trends. 
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